In [1]:
%pylab inline
import pandas as pd
from scipy import stats
import statsmodels.api as sm
import matplotlib.pyplot as plt
import warnings
from itertools import product
def invboxcox(y,lmbda):
if lmbda == 0:
return(np.exp(y))
else:
return(np.exp(np.log(lmbda*y+1)/lmbda))
In [2]:
data = pd.read_csv('WAG_C_M.csv',';', index_col=['month'], parse_dates=['month'], dayfirst=True)
data.columns = ['Salary']
plt.figure(figsize(15,7))
data.Salary.plot()
plt.ylabel('Salary')
pylab.show()
У ряда есть явный тренд и у него разная дисперсия. Проведем формальную проверку + STL разложение
In [3]:
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.Salary).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.Salary)[1])
Сделаем преобразование Бокса-Кокса для стабилизации дисперсии:
In [4]:
data['Salary_box'], lmbda = stats.boxcox(data.Salary)
plt.figure(figsize(15,7))
data.Salary_box.plot()
plt.ylabel(u'Transformed Salary')
print("Оптимальный параметр преобразования Бокса-Кокса: %f" % lmbda)
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.Salary_box)[1])
plt.show()
У нас есть выраженная сезонность и тренд. Пробуем сделать диффиренцирование на сезонность
In [5]:
data['Salary_box_diff'] = data.Salary_box - data.Salary_box.shift(12)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.Salary_box_diff[12:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.Salary_box_diff[12:])[1])
plt.show()
Хотя гипотеза нестационарности отвергается критерием Дики-Фуллера, тренд все еще есть. Нужно делать дифференцирование.
In [6]:
data['Salary_box_diff2'] = data.Salary_box_diff - data.Salary_box_diff.shift(1)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.Salary_box_diff2[13:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.Salary_box_diff2[13:])[1])
plt.show()
Можно проверить как выгледит ряд при другом порядке диффиринцирования.
In [7]:
data['Salary_box_diff_invert'] = data.Salary_box - data.Salary_box.shift(1)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.Salary_box_diff_invert[1:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.Salary_box_diff_invert[1:])[1])
plt.show()
In [8]:
data['Salary_box_diff_invert2'] = data.Salary_box_diff_invert - data.Salary_box_diff_invert.shift(12)
plt.figure(figsize(15,10))
sm.tsa.seasonal_decompose(data.Salary_box_diff_invert2[13:]).plot()
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(data.Salary_box_diff_invert2[13:])[1])
plt.show()
Оба варианта приводят к стационарным рядам и можно начинать подбирать начальные параметры.
In [9]:
plt.figure(figsize(15,8))
ax = plt.subplot(211)
sm.graphics.tsa.plot_acf(data.Salary_box_diff2[13:].values.squeeze(), lags=48, ax=ax)
pylab.show()
ax = plt.subplot(212)
sm.graphics.tsa.plot_pacf(data.Salary_box_diff2[13:].values.squeeze(), lags=48, ax=ax)
pylab.show()
Начальные приближения: Q=0(последний значинимы сезонный лаг / 12), q=6(последни значиммый не сезонный лаг, но не сильно большой), по частичной автокоррекции подбираем P=1(последний значинимы сезонный лаг / 12), p=6(последни значиммый не сезонный лаг, но не сильно большой)
In [10]:
ps = range(0, 7)
d=1
qs = range(0,7)
Ps = range(0, 2)
D=1
Qs = range(0,2) # Т.к. мы будем подбирать лучшую модель по гиперпараметрам, позволим немного расширить варианты подбора задав Q=1
parameters = product(ps, qs, Ps, Qs)
parameters_list = list(parameters)
len(parameters_list)
Out[10]:
In [11]:
%%time
results = []
best_aic = float("inf")
warnings.filterwarnings('ignore')
for param in parameters_list:
#try except нужен, потому что на некоторых наборах параметров модель не обучается
try:
model=sm.tsa.statespace.SARIMAX(data.Salary_box, order=(param[0], d, param[1]),
seasonal_order=(param[2], D, param[3], 12)).fit(disp=-1)
#выводим параметры, на которых модель не обучается и переходим к следующему набору
except ValueError:
print('wrong parameters:', param)
continue
except LinAlgError:
print 'wrong parameters (LinAlgError):', param
continue
aic = model.aic
#сохраняем лучшую модель, aic, параметры
if aic < best_aic:
best_model = model
best_aic = aic
best_param = param
results.append([param, model.aic])
warnings.filterwarnings('default')
In [12]:
result_table = pd.DataFrame(results)
result_table.columns = ['parameters', 'aic']
print(result_table.sort_values(by = 'aic', ascending=True).head(10))
Минимальный AIC мы получаем для модели с параметрами (5,5,1,0). Модели с более маленькими коофициентами имеют существенную разницу в AIC, как лучшую модель берем именно ее. Посмотрим на параметры лучшей модели
In [13]:
print(best_model.summary())
In [14]:
plt.figure(figsize(15,8))
plt.subplot(211)
best_model.resid[13:].plot()
plt.ylabel(u'Residuals')
ax = plt.subplot(212)
sm.graphics.tsa.plot_acf(best_model.resid[13:].values.squeeze(), lags=48, ax=ax)
print("Критерий Стьюдента: p=%f" % stats.ttest_1samp(best_model.resid[13:], 0)[1])
print("Критерий Дики-Фуллера: p=%f" % sm.tsa.stattools.adfuller(best_model.resid[13:])[1])
Остатки несмещены (подтверждается критерием Стьюдента) стационарны (подтверждается критерием Дики-Фуллера и визуально), неавтокоррелированы (подтверждается критерием Льюнга-Бокса и коррелограммой). Посмотрим, насколько хорошо модель описывает данные:
In [15]:
data['model'] = invboxcox(best_model.fittedvalues, lmbda)
plt.figure(figsize(15,7))
data.Salary.plot()
data.model[13:].plot(color='r')
plt.ylabel('Salary')
pylab.show()
Визуально модель описывает данные хорошо.
In [16]:
data2 = data[['Salary']]
date_list = [datetime.datetime.strptime("2016-08-01", "%Y-%m-%d") + relativedelta(months=x) for x in range(0,36)]
future = pd.DataFrame(index=date_list, columns= data2.columns)
data2 = pd.concat([data2, future])
data2['forecast'] = invboxcox(best_model.predict(start=284, end=284+36), lmbda)
plt.figure(figsize(15,7))
data.Salary.plot()
data2.forecast.plot(color='r')
plt.ylabel('Salary')
pylab.show()
Визуально модель делает достаточно хороший прогноз.
In [ ]:
In [ ]: